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' ^ ■ Abstract 

This technical report yields detailed calculations of the paper [JJ which have been however 
••O | automated since (see [2]). It deals with the stability analysis of various finite difference schemes 

for Maxwell-Debye and Maxwell-Lorentz equations. This work gives a systematic and rigorous 
| continuation to Petropoulos previous work [B]. 

r^* ! 1 Introduction 
-i— > 

c3 . 

We address the stability study of finite difference schemes for Maxwell-Debye and Maxwell-Lorentz 
models. To this aim we selected the same schemes as those already studied by Petropoulos [6j, who 
after having correctly denned characteristic polynomials associated to each scheme, merely computed 
its roots with a numerical algorithm. This implies having to specify values for the physical parameters 
ON I which occur in the models as well for the time and space steps chosen for the discretization. The 
analysis has therefore to be carried out anew for each new material or discretization. We perform here 
a von Neumann analysis on the characteristic polynomials in their literal form, which yields once and 
\ for all stability conditions which are valid for all materials. 

o 

q ■ 1.1 Maxwell-Debye and Maxwell-Lorentz Models 



- 

X 



Le us consider Maxwell equations without magnetisation 



(Faraday) 9 t B(i,x) = - curl E(t,x), 

(Ampere) d t D(i,x) = — curlB(t,x), ^ 

Mo 

where x 6 H. . This system is closed by the constitutive law of the material 

D(t,x) = e £ooE(t,x)+£o f E(t-r,x)x(r)cfr, (2) 

JO 

where £oo is the relative permittivity at the infinite frequency and x the linear susceptibility. If we 
discretize the integral equation ([2]), we obtain what is called a recursive scheme (see e.g. [5], [TO]). 
We can also differentiate Eq. (|2|) to obtain a time-differential equation for D which depends on the 
specific form of x- F° r a Debye medium, this differential equation reads 

t r dtD + D = t r e £oo^E + e £sE, (3) 

where t r is the relaxation time and e s the static relative permittivity. We can derive an equivalent 
form dealing with the polarisation polarisation P(i,x) = D(t,x) — eo £ ooE(t,x), namely 

t r 5 t P + P = e (es-£oo)E. (4) 
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For a Lorentz medium with one resonant frequency u)\, we have similarly 

d?T> + udtD + = e £oo^ 2 E + e £oo^E + eo^fE, (5) 

where v is a damping coefficient and 

d?P + udtP + ujfP = e (e s - £oo)a;?E. (6) 

Denoting by J the time derivative of P, Maxwell system ([1]) can be cast as 

&B(t,x) = -curlE(i,x) 
1 

1.2 Yee Scheme 



1 (7) 
£o£oo<%E(i, x) = — curl B(i, x) — J(t,x). v ; 



To discretize Maxwell equations in a passive medium (J = 0), we use Yee scheme [8], which consists in 
stagg ering space and time discretization grids for the different fields. We denote by Cqq — 1/ yj £0^-00/^0 
the light speed at infinite frequency. If the space step 5x is the same in all directions and 5t is the time 
step, the CFL condition is c^St/Sx < 1 in space dimension N = 1 and Coodt/dx < l/v2 for N = 2 or 3. 
In dimension 1, we can for example only consider fields E = E x et B = B y which discrete equivalents 

are E? ~ E(n5t, j5x) (with similar notations for D = D x ) and B n+ f ~ B((n + h)5t, (j + h)5x). Yee 

3+2 

scheme for the initial Maxwell system (pQ) in variables E, B and D therefore reads 



1 ,_„.i „„, 1 ._n+i _n+i. W 



5t y 3 3 /x <5x 3-5 

In the same way, for Maxwell system ([7]) in variables E, B and J, we have the Yee discretization 

1 {B n+h _ ^n-i = 1 ■ EJ), 

St J+2 J+2 dx 3 3 (9) 

^{EZ* 1 - EP) = -J-(B n +? - B n+ ?) - J" + * ■ 

dt 3 J HqOX 1+2 3~2 3 

For the matter equations, we address "direct integration" schemes which discretize the differential 
equations ©-© (see @], [3], [9]). 

Before describing and analysing the schemes one by one, we give below the principle of the von 
Neumann analysis which allows us to study their stability. 

2 Principles of the von Neumann Analysis 
2.1 Schur and von Neumann polynomials 

We define two families of polynomials: Schur and simple von Neumann polynomials. 
Definition 1 A polynomial is a Schur polynomial if all its roots r satisfy \r\ < 1. 

Definition 2 A polynomial is a simple von Neumann polynomial if all its roots r belong to the unit 
disk (\r\ < 1) and all the roots of modulus 1 are simple roots. 
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It may be difficult to localise roots of a polynomial with complicated coefficients. On the other 
hand, we can turn this difficult problem into the solving of many simpler small problems. To this aim, 
we construct a polynomial series with strictly decreasing degree. To a polynomial <p defined by 

4>(z) = c + c\z H h c p z p , 

where Co, c\ . . . , c p £ C and c p 7^ 0, we associate its conjugate polynomial 4>* which reads 

4>*{z) = c* + c* v _ x z + ■■■ + c* z p . 

Given a polynomial 4>o, we can define a series of polynomials by recursion 

, / \ _ <t)*m{Q)<Pm{z) ~ <pm(0)<P* m (z) 
<Pm+l\Z) — • 

Z 

This series is finite since it is clearly strictly degree decreasing: deg</> m +i < deg</> m , if 4> m ^ 0. Besides, 
we have the following two theorems at our disposal. 

Theorem 1 A polynomial (fr m is a Schur polynomial of exact degree d if and only if 4> m +i is a Schur 
polynomial of exact degree d — 1 and |0 m (O)| < |^(0)|. 

Theorem 2 A polynomial <j) m is a simple von Neumann polynomial if and only if 
4>m+i is a simple von Neumann simple polynomial and \4> m (0)\ < |^^(0)|, 

or 

4>m+i is identically zero and <f>' m is a Schur polynomial. 

To localise roots of rf>o in the unit disk or not, we only have to check conditions at each step m (non 
zero leading coefficient, |(/> m (0)| < |^(0)|, . . . ) until we obtain a negative answer or a polynomial of 
degree 1. 

The proofs of the above results are based on Rouche theorem and are given in [7] . 



2.2 Stability Analysis 

The models we consider are linear. They can therefore be analysed in the frequency domain. Hence 
we assume that the scheme deals with a variable with space dependency in the form 

^ = C/"exp(^-j), 

where £ et j S R , iV = 1, 2, 3. Let G be the matrix such that U n+1 = GU n and we assume it does not 
depend on time, nor on 5x and 5t separately but only on the ratio 5x/5t. Let <j>o be the characteristic 
polynomial G, then we have the following sufficient stability condition. 

Theorem 3 A sufficient stability condition is that 4>o is a simple von Neumann polynomial. 

This condition is not a necessary one. The stability is linked to the fact that U n = G n U° and 
corresponds to the boundedness of the iterates G n of the matrix G. The case of multiple unit modulus 
roots can give rise to iterates of G which are bounded (e.g. for the identity matrix) or not. For 
example 

^ ^ is bounded, and [ \ \ ) = ( \ ^jis not bounded. 



1 / V 01 / ' V 01 / V 1 , 

This case occurs for the schemes we are dealing with and have to be treated separately, without the 
help of the von Neumann analysis, which handles characteristic polynomials and not the matrices they 
stem from, which induces a loss of information. 
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3 Debye Type Media 



For Debye type media, we study two schemes. The first one is due to Joseph et al. [3] and consists 
in coupling Maxwell equations in variables E, B and D with the Debye model linking E and D. The 
second is due to Young [9] and couples Maxwell equations in variables E, B and J with the Debye 
model linking E, P and J. 



3.1 Joseph et al. Model 
3.1.1 Model Setting 

Maxwell system ([8]) is closed by a discretization of the Debye model 



E 



St 



e™ +1 + m 

+ e e s ^— 3 - 



D 



n+l 



tr 



St 



System 



TUl) deals with the variable 



namely 



(10) 



n\t 



and reads 



3 + i 



B 



i+2 J+5 



Sx 



\pj+i z-j 



OX 3 + 2 



£ 



n+l 



We see that this formulation contains dimensionless parameters: 

CFL constant, 



,n+l 



A = CooSt/Sx 
5 = St/2t T 



-oo 



normalised time step, 
normalised static permittivity. 



We write this system into the explicit form 



n+l 



3+h 



A(^ +1 -e?), 



v 

j 

(l + Se> s )£] +1 



(1 - Se' a )£« + (1 + S)X 2 (£j L +1 - 2£« + q_ x ) 



which yields the amplification matrix 



G 



( 1 -A(e* - 1) 

(l+5)X(l-e- i i) (l-feQ+(l+<5)A 2 (e^-2+e-^) 

A 2 (e^ -2 + e~^) 





25 
l+5e' B 

1 



We set a 
reads 



A(e 



\ "A(l 

1) and q = \a\ 2 = -A 2 (e^ - 2 + e - *) = 4A 2 sin 2 (£/2). With these notations G 



G 



( 1 

(ljjK (l-fe£)-(l+% 
l+5e' s l+5e' s 

-q 




1 
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3.1.2 Computation of the Characteristic Polynomial 

The characteristic polynomial of G is equal to 



P{Z) 



1 



l + Se' s 
1 



Z-l a 

-(l + 8)a* (1 + 8e' s )Z - (1 - Ss' s ) + (1 + 6)q -25 
-a* q Z-l 



l + 5e: 
1 



7 (S-l) 



1 + fe'jZ - (1 - 5e' s ) + (1 + <J)g -25 
g Z-l 



1 + 5 -5 
1 Z-l 



j + fe 7((^ + fesl^ 2 " P " (1 + <W + [1 " 54 - (1 - 5)q]} 

+q {[l + 6]Z-[l-5]} 
The characteristic polynomial is proportional to 

M z ) = [i + fe' s ]^ 3 - [3 + 54 - (i + %]^ 2 + [3 - 54 - (i - - [i - K}- 

3.1.3 Von Neumann Analysis 

From the polynomial 0q, we perform the recursive construction of the above-mentioned series of 
polynomials. We therefore define 

4>l(Z) = [1 + 5e' s ] - [3 + fe' s - (1 + S)q]Z + [3 - <fe' s - (1 - 5)q}Z 2 - [1 - fe^Z 3 . 

The condition |</>o(0)| < |0g(O)| is valid. We define by recursion 

= 25{2e' s Z 2 - [44 - (4 + l)q]Z + [24 - (4 - l)q]}, 
<t>{(Z) = 25{24-[44-(4 + l) (? ]Z+[24-(4-l)g]Z 2 }. 

Since e s > e^, we have e' s > 1 and the quantity e' s — 1 is nonnegative. If q = or 4 = lj we have 
exactly |0i(O)| = |0*(O)|, and these specific cases have to be treated separately (see below). In the 
opposite case, condition |<^>i(0)| < |0*(O)| reverts to (e' s — l)q < Ae' s . It is reasonable to assume we will 
not obtain a better result than with the raw Yee scheme (A < 1) and therefore q € [0, 4]. In that case, 
and provided q / 0, we do have |^i(0)| < |<^(0)|. Moreover the degree of polynomial <f>\ is 2. Last 

MZ) = ^Uo)Mz) - M0)4>UZ)) 

= 4<5 2 (4 - l)g[(44 - (4 - l)q)Z ~ (44 - (4 + 1)?)] • 
Always in the case when e' s > 1 and q E]0, 4], the leading coefficient 

44-(4-i)( ? = 4(4- (? ) + ^o, 

and is degree of fa is 1. The root of </>2 is 

v 44 - (4 + l)q 



4e' 



1) 



are Schur 

1 are simple 



The modulus of this root is strictly lower than 1 if q ^ 4 and therefore 4>2 and hence 
polynomials thanks to Theorem [TJ If q = 4, the root of 4>2 is —1 and </>2 and hence 
von Neumann polynomials thanks to Theorem [2j In both cases, we obtain the stability with the only 
assumption A < 1, provided we treat the above-mentioned special cases. 
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3.1.4 Case q = 



The case when q = corresponds to the characteristic polynomial 

MZ) = (Z - 1) 2 ([1 + 5e' s ]Z - [1 - 5e' s ]) 

which is not a simple von Neumann one. We shall therefore study the amplification matrix directly, 
which is then simply 

/ 1 \ 

1 - 8e' a 25 





G 



l + 5e' s 




1 + Se' s 
1 



/ 



We clearly see that the eigenvectors corresponding to the eigenvalue 1 are in two stable eigensubspaces. 
The other eigenvalue has a modulus strictly lower than 1. Iterates of this matrix are therefore bounded. 
This conclusion is valid for all e' > 1. 



3.1.5 Case e' s = 1 

The case when e' s = 1 gives rise to a breaking of condition |</>i(0)| < |0*(O)|. We shall therefore study 
directly the nature of <f>\ without carrying recursion over. We have 

MZ) = 45{Z 2 -[2-q]Z + l}, 

which determinant is q(q — 4) and is therefore negative if q G]0,4[. The roots of 4>\ are therefore 
complex conjugate, distinct and their modulus is 1 (their product is equal to 1). The polynomial 4>i 
is therefore a simple von Neumann one, and 4>q also. 



3.1.6 Case q = 4 

The last case we have to treat is e' s = 1 and q = 4, where — 1 is a double root of 4>\ . It is also a double 
root of 0o which reads 

MZ) = {Z + lf([l + 5}Z-[l-5}). 
Hence we study directly the amplification matrix which reads simply 




There is no trivial splitting in two distinct eigensubspaces. We compute the eigenvectors associated 
to the eigenvalue —1. To this aim we solve 

/ 2 -a 0\ / 2 -a \ 

(G + Id)F= U* if,-? S 7 = 0^ 1 -1 ]V = 
\a* -q 2 J V a* -q 2 J 

and we only find one eigendirection, that of V = (a, 2, 2)'. A minimal two-dimensional eigensubspace 
is therefore associated to the eigenvalue —1 and iterates G n are linearly increasing with n. Hence we 
conclude to instability when q = 4 and e' s = 1. 
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3.1.7 Synthesis for the Debye Joseph et al. Model 

The scheme (j8])- (|10p for the one-dimensional Maxwell-Debye equation is stable with the condition 

St < 8x/ Cqo if £s > £00 and 5t < 5x/ if e s = e^. 

We have already seen in this first example different types of arguments to conclude to stability: the 
generic case {q G]0,4[ and e s > 1) gives rise to a Schur polynomial via Theorem [IJ the cases q g]0, 4[ 
and e' s = 1 or q = 4 and e' s > 1 to a simple von Neumann polynomial via Theorem [2] and last, the 
case q = to a (not simple) von Neumann polynomial, but with a double eigenvalue that operates on 
two stable and distinct eigensubspaces. We have also encountered an instable case when e' s = 1 and 
q = 4 which nevertheless corresponds to a (non simple) von Neumann polynomial. 



3.2 Young Model 
3.2.1 Model Setting 

Maxwell system ([9]) is closed by two discretizations of Debye equation namely 



St 

and 



p 1 2 p 2 p 1 2 _|_ p' 

TT ~ = ~- + Z0(e s ~ Eoo)E?, (11) 



t r 2 = -Pj 2 + E {e s - £oo) — — - -. 



:i2) 



71 H™ — 

Although we make use of J- 2 in the description of the scheme, this is not a genuine variable and the 
system (f^- (fTT]l - (fT2l) deals with the variable 



and reads 



J T 2 J T 2 



In this system, apart from the notations A, 5 which we have already defined, we have introduced the 
dimensionless parameter a = e' s — 1 which is , as we already mentioned, a non negative parameter. 
We rewrite this system in the explicit form 



[l + 5a)£^ = (i-S a )£?-\(B':j-B i :j) + \<(£? +1 -2£2 + £^ 1 ) 



?n + i , | ,;„,o« \(K n 2 _ Q n 2 N 

3+5 J' 2 

+2<5^V^ + 



(1 + 5)7^ 2 = (1 - S)V- 2 + 25a£f, 
from which stems the amplification matrix 

/ 1 -a 



G 



{l+8)(l-8a)+A8 2 a-(l+8)q 1-8 28 

l+8a (1+5)(1+Sa) l+<5 l+8a 

28a 1-8 



1+5 1+8 
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3.2.2 Computation of the Characteristic Polynomial 

The characteristic polynomial G is 



P(Z) 



l+Sa 




1 a 

(1+8) (l-8a)+A8 2 a- (l+8)q 
(l+S)(l+8a) 



z- 



28a 

T+3 





1-8 28 
1+8 l+8a 

z- 1 - 5 



To reduce computations, we set Y = Z — 1, which yields 
{\ + 5) 2 {l + 5a)P{Z) 



Y a 

-(1 + <5)ct* (l + 6)(l + 8a)Y + 26a(l-S) + (l + S)q -25(1-5) 
-25a (1 + 5)Y + 25 

Y a 
-(l + 5)a* (1 + 5)(1 + 5a)Y + (1 + 5)q (1 + <5)(1 - 6)Y 

-25a (1 + 5)Y + 25 



We see that (1 + 5) is a factor in both sides and therefore 
(1 + 5)(l + 5a)P(Z) = Y 



(l + 5a)Y + q (1-5)Y 
-25a (1 + 5)Y + 25 



+ q 



1 -25(1-5) 
(l + 5)Y + 25 



= Y {[(1 + + 5a)]y 2 + [25(1 + q) + (1 + %]Y + [25g]} 

+g(l + 5){[l + <5]Y + [25]} 
= [(1 + <5)(1 + <fo)]Y 3 + [2(5(1 + a) + (1 + %]Y 2 

+ [(l + 35)q]Y +[25q}. 

The characteristic polynomial is proportional to 

MZ) = [{l + 5a){l + 5)]Z 3 -[3 + 5 + 5a + 35 2 a- {l + 5)q]Z 2 
+[3 - 5 - 5a + 35 2 a - (1 - 5)q]Z - [(1 - 5a)(l - 5)]. 

3.2.3 Von Neumann Analysis 

Condition |</>o(0)| < |<$)(0)| is valid without any assumption. We define by recursion 

M z ) = 25{[2(l+a){l + 5 2 a)]Z 2 - + a){l + 5 2 a) - {2 + a + 5 2 a)q]Z 
+ [2(1 + a) (1 + 5 2 a) - a(l - 5 2 )q]}. 

The case when 5 2 > 1 does not allow to fulfil the condition |</>i(0)| < |0*(O)|. We will assume therefore 
for the von Neumann analysis that 5 < 1, which bounds the time step with respect to the time delay 
t T . This is reasonable from the point of view of modelling: we cannot approximate the delay equation 
with too large a time step. Such an assumption was however not necessary for the Joseph et al. 
scheme. 

The equality case |0i(O)| = |0i(O)| is obtained when q = 0, a = (i.e. e' s = 1) or 5 = 1. These cases 
shall be treated separately again. 

If a > 0, q > and 5 < 1, then \(f>i(0)\ < \</>t(0)\ is equivalent to a(l-5 2 )q < 4(1 + a)(l + 5 2 a), which 
is clearly true if q g]0,4]. Besides, the degree of polynomial 4>i is 2. 
In the general case (a > 0, q G]0,4] and 5 < 1), we then compute 02 

<fo(Z) = A5 2 a{l-5 2 )q{[A(l + a)(l + 5 2 a)-a{l-5 2 )q]Z 

-[4(1 + a)(l + 5 2 a) - (2 + a + 5 2 a)q}}] . 
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We split the study according to the sign of 2 (0) (02 (0) is clearly always positive) and in both cases 
|02(O)| < |(/>2(0)|, for q e]0,4]. The root of 02 therefore belongs to the interval ] — 1,1[ and 0o is a 
Schur polynomial. Hence we obtain the stability with the assumptions A < 1 and 5 < 1, provided we 
treat the above-mentioned specific cases. 



3.2.4 Case q = 

The case when q = corresponds to the characteristic polynomial 

Mz) = (z-i) (z- {l + s){1 + Sa) ) 

and is not a simple von Neumann one. The amplification matrix reads 



G 



( 1 







(l+<5)(l-<5a)+4<5 2 q 

(1+S)(l+Sa) 
25a 
1+5 





1-5 25 
1+5 l+5a 

1-5 

1+5 



We clearly see that the eigenvectors corresponding to the eigenvalue 1 are in two stable eigensubspaces. 
Iterates of this matrix are therefore bounded. This conclusion is once more valid in the limit cases 
8 = 1 and a = 0. 



3.2.5 Case e' s = 1 

We notice that 0o is the same as that for the Joseph et al. model for e' s = 1. Polynomial 0o is therefore 
a simple von Neumann polynomial for 7^ 4 (see above). The value of 5 does not play any role here. 
This corresponds to different amplification matrices, operating on different sets of variables, the link 
between both formulations being not straightforward. Case q = 4 has therefore to be treated anew. 



3.2.6 Case q = 4 

If q = 4, only the case when a = has not been treated by the general study and —1 is once more a 
double eigenvalue 

-a 

rl-<5 



(G + Id) V =| a* 2-q 2<5^§ \ V = 0^\a*2-q | V = 0, 

T+5 



2 




and the only eigendirection is that of V = (a, 2, 0)', which gives rise to linearly increasing iterates G r ' 
and to instabilities. The fact that 5 = 1 or not does not play any role in this argument. 



3.2.7 Case 5 = 1 

There remains to study the case S = 1 for q G]0, 4] and a > 0. Then Z = is a trivial root of 0o, 
which simply reads 

MZ) = 2(1 + a)Z{Z 2 - [2 - -^—\Z + 1}. 

I + a 

The discriminant of the second order factor is A = fjj^p (<7 — 4(1 + a)) < 0. We therefore have 
two distinct complex conjugate eigenvalues of modulus 1. Polynomial 0o is a simple von Neumann 
polynomial. 
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3.2.8 Synthesis for the Debye Young Model 

The scheme ([9|)- (jlip - (|12p for the one-dimensional Maxwell-Debye equation is stable with the condition 

St < min(fe/c 00 , 2t T ) if e s > and St < dx/coo if e s = £oo- 

If e s > £oo, the stability condition is more restrictive for the Young scheme than for Joseph et al. 
scheme. The obtained bound is also related to the good approximation of Debye equation. 



4 Anharmonic Lorentz Type Media 

For Lorentz type media, we study three schemes. The first one is due to Joseph et al. [3] and consists 
in coupling Maxwell equations in the variables E, B and D with the Lorentz model linking E and 
D. The second and third are due to Kashiwa et al. [4] and Young [9] respectively and both couple 
Maxwell equations in the variables E, B and J with the Lorentz model linking E, P and J. They 
differ in the choice of the time discretization of J. 



We restrict here to the anharmonic case for which the damping v is non-zero. The harmonic case 
{y = 0) is treated with the same schemes but the analysis happens to be much more technical. To 
keep proofs readable in the general case we postpone the harmonic case to the next section. 



4.1 Joseph et al. Model 
4.1.1 Model Setting 

Maxwell system ([8]) is closed by a discretization of the Lorentz equation (|5|), namely 
E n+1 — 2E n + E n ~ x E n+l — E n ~ l E n+1 + E n ~ x 

3 3 3 i , , — . 3 3 i — — . , ,2 3 3 



£q£oc r- vsq£oo ^ e o £ s U) i 



5t 2 u * 2St 
D n+i _ 2D™ + D n ~ x D" +1 - B n ~ x „ D" +1 + D"- 1 
+ v— — ^r-^ Voj{— — . 



(13) 



5t 2 25t 

The explicit version of system ([5|)- (fT3"]) does not use explicitly the variable D™~ . Indeed we can use 
the system deals with the variable 



the explicit formula to compute D™ + — D™ and the implicit one to compute — and therefore 



and reads 



i 



n \ 



(1 + 5 + u;e' s )q +l = 2£] + (1 + S + io)\\q +l - 2£? + £"_ x ) -{1-5 + coe'^q- 1 



3 

l „ l 



-25X(B I -B 1 2 ) + 2wP™. 

In this system, apart from the already used notations A and e' s , we have denoted 

5 = 5ti//2 normalised time step, 

uj = u)\5t 2 /2 square of the normalised frequency. 
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The amplification matrix of the system is 





1 


—a 








\ 




25a* 


2-q(l+5+ui) 


l-S+we' s 










l+S+uie' s 


l+S+we' e 











1 










V 


a* 


-q 





1 


/ 



G 



4.1.2 Computation of the Characteristic Polynomial 

The characteristic polynomial of G is equal to 



P{Z) 



Z - 1 

2&o* 



-a* 



Z - 



a 

2-q(l+<S+cj) l-8+ue' s 





2k; 
1+<5+c< 





Therefore 

(l + 5 + ue> s )P{Z) 

Z - 1 



-2(5a* (l + <5 + u;4)Z-2 + (l + <5 + u;)g l-fl + we^ 
-1 



-a 



+ Z 



Z - 1 






z - 










-2w 


Z 








Z-l 




a 



-25a* (1 + <5 + uje' s )Z - 2 + (1 + 6 + u)q -2u 
-a* q Z-l 



Z-l 

-25a* 1 - 5 + ue' s -2uj 
-a* Z-l 

= {Z - l) 2 [l - 5 + ue' s ] 

+Z{Z - 1){[(1 + (5 + oje'jZ - 2 + (1 + 5 + w)q] [Z - 1] + 2wg} + g Z[25(Z - 1) + 2w]. 

The characteristic polynomial is proportional to 

(j)o(Z) = [I + 5 + u;e' s }Z 4 - [4 + 25 + 2u;e' s - (I + 5 + cu)q]Z 3 + [6 + 2^4 - 2q]Z 2 
- [4 - 25 + 2ue' s -{1-5 + u)q}Z + [1 - 5 + Loe' s }. 







4.1.3 Von Neumann Analysis 

We successively compute 

MZ) = 2<5{[2(1 + cje' s )]Z 3 - [6 + 4cje' s - (2 + u{l + e' s ))q]Z 2 

+ [6 + 2coe' s - 2q\Z - [2 + u{e' s - l)q]}, 
MZ) = 45 2 io{[4e' s (2 + Loe' s )-4(e' s -l)q-u{e' s -lfq 2 ]Z 2 

- [84(2 + ue' s ) - 4((4 - 1) - 4(2 + ue' s ))q + 2(4 - l)q 2 ]Z 
+ [44(2 + ue' s ) - 4(4 - 1)(2 + u;e' s )q + (2 + u{l + 4))(4 - l)g 2 ]}, 
fo(Z) = U5W(e' s -l)(l + u;e' s )q(2-q)x 

x {[44(2 + ue' s ) - (4 - 1)(6 + 2ue' s ))q + (4 - 1)(1 + ^)g 2 ]^ 

- [44(2 + coe' s ) - 2((4 - 1) + 4(2 + ue' s ))q + (e' s - l)q 2 ]}. 

The root of fa is 

z= (2- g )(2£ / s (2 + g;4)-(4-l) g ) 

(2 - <z)(24(2 + ue' s ) - (e> - l)g) + 2(2 + uje' s )q + (e' s - 1)(1 + w)g 2 ' 
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form on which we easily see that Z remains of modulus < 1 if g g]0, 2[, which corresponds to the 
condition for a multi-dimensional Yee scheme (A < l/\/2). We notice from now on that we shall treat 
cases q = and q = 2 apart because 03 = 0. In the general case, we have to check the intermediate 
properties. First, the degree of polynomials 0i and 02 is 3 and 2 respectively. The degree of polynomial 
03 is 1 provided e' s > 1. We shall treat the case e' s = 1 apart. In the general case (q/0 and e' s > 1), 
there remains to check the estimates between 

0o(O) = 1-<5 + w4, 

$(0) = l + 5 + u,4, 

01 (0) = 2<5[-2-o,(4-l)< ? ], 

0i(O) = 25[2(lW B )], 

0a(O) = 4<5 2 u,[44(2 + o<) - 4(4 - 1)(2 + we'Jg + (2 + co(l + 4))(4 - 1)^] , 

02(0) = 45 2 co [44(2 + 0,4)- 4(4 - i) g - ^ - 1)V] • 

It is clear that for e]0, 2[, we have |0o(O)| < |0o(O)| and |0i(O)| < |0i(O)|. A simple calculation shows 
that 

02(0) = 4(5 2 o;[(2-g) 2 (2(l+a;)(4-l)+a;(4-l) 2 ) + 8 + 4a; + 4a;(4-l)g], 
05(0) = 4^[(2- (7 )(4(4-l)+ W (4-l) 2 g) + 8 + 4^ + 8^(4-1)], 

form on which we readily see that both quantities are positive. Besides 

02(0) - 02(0) = 85 2 oj(1 + u,4)(4 - l)q(2 -q)>0. 

We therefore checked all the assumptions. In the general case, 0o is a Schur polynomial. 



4.1.4 Case q = 

The case q = gives anew rise to a separate study. We have 

O (Z) = (Z-l) 2 [(Z-l) 2 + <5(Z 2 -l)+o,4(Z 2 + l) 

The corresponding amplification matrix is 

/ 1 





2 





1-5 + oje' 





2o, 



\ 



1 + 5 + ue' s 
1 

V o o 



1 + 5 + Lue' s l + 5 + Luei 










1 



/ 



Once more, 1 is a double root but in two distinct eigensubspaces. There remains to check that the 
other factor of the polynomial, namely 

ifo(Z) = [1 + 5 + coe' s }Z 2 - 2Z + [1 - 5 + ue' s ], 

is a Schur (or a simple von Neumann) one. We do have | -0o(O)| < (V'o(O)l and we compute 

MZ)=^{[l+ue' s ]Z-l}. 

The modulus of both remaining eigenvalues is strictly less that 1 and iterates of the amplification 
matrix are bounded. This holds whatever the value of e' 
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4.1.5 Case q = 2 

In the specific case q = 2, 0o reads 

</> {Z) = [1 + 6 + ujs' s }Z a - [2 + 2co(e' s - 1)]Z 3 + [2 + 2uje' s ]Z 2 - [2 + 2uo{e' s - 1)]Z + [1-5 + Loe' s ] 
= (Z 2 + 1){[1 + 5 + ue' s ]Z 2 -[2 + 2co(e' s - 1)]Z + [1 - 5 + ue' B ]}, 

which has +i as simple roots. This is therefore a good candidate to be a simple von Neumann 
polynomial. The remains to study the other factor of the polynomial 



which is a Schur polynomial for all e' s . Polynomial ipo is therefore a Schur polynomial and 4>o is a 
simple von Neumann polynomial. 

4.1.6 Case e' s = 1 

If £g = 1, the polynomial 03 is identically zero and 



The root of 2 does have a < 1 modulus if q e]0, 2]. The polynomial 0o is therefore a simple von 
Neumann polynomial. 

4.1.7 Synthesis for the Lorentz Joseph et al. Model 

The scheme (f8|)- (|13p for the one-dimensional anharmonic Maxwell-Lorentz equations is stable with 
the condition 



4.2 Kashiwa et al. Model 
4.2.1 Model Setting 

A modified version of Maxwell system ([9]) is closed by a discretization of Lorentz equation ([6]) , namely 



MZ) = [1 + 6 + ue' s }Z 2 -[2 + 2uj(e' s - 1)]Z + [1 - S + us' s ] 
which has not +i as roots. We notice that |V>o(0)| < IV'oWI arm compute 

MZ) = 4<5{[1 Wjz - [1 + ^(4 - 1)]}, 



Mz) 
Uz) 



165 2 w(2 + oj){Z 2 - [2 - q]Z + 1} 
16S 2 u{2+u){2Z - [2-q]}. 



5t < 5x/V2> 





- v -(jf l + j?) + 



U)\(e s - £oo)g 

2 



(ej +i + e?) - + if). 



(14) 



The system (I14p deals with the variable 



2 cn -pn qn\t 
\ ' °3 ' ' 3 ' °3 > 



and reads 
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[l + 5+ l -ue' s }£^ 1 = [i + 5- l -u{e' s -2)]£^ + \\l + 5+ l -u)(£^ 1 -2£^ + £^ 1 ) 

_ A (i + s + \^)(B r -ii -q:p+tsp?-j?, 
[i + s + \ue' s \v] +i = [i+5+^(4-2)]p™-^4-i)(^;;|-fi;_i) 

+u(e' s - 1)£* + \\Me's ~ l)(^+i - 2£? + fJU) + JJ\ 

[l + tf + ^jjj* 1 = [i-i-^-Mei-i)(B^J-r|) 

from which stems the amplification matrix 



G = 



I 


1 




— 0" 












CT *(D-i W (4-l)) 


(l-q)D 


-(2-,)|w(ei-l) 


Ul 


-1 






r> 

<T^w(ei-l) 


(2- 


D 

«)J«(ei-l) 


D 
D-u) 


D 
1 




V 


L> 
D 


(2- 


-?M4-i) 


D 
-2w 
D 


D 

2-D 
D 


J 



where, along with earlier notations, D = 1 + 5 + \ive' s 



4.2.2 Computation of the Characteristic Polynomial 

The characteristic polynomial of G is 



P(Z) 



Z-l 

-a*(D-±u,(e' s - 



1)) 



i) 



D 
D 



hence setting X = D(Z — 1) 



Z- 



(l-q)D-(2~q)^(e' s -l) 
D 

(2-9)iw(e;-l) 



i) 





_ 0£ 
D 

D-ui 



D 



2uj_ 
D 





J_ 

D 

-1 
D 

rj 2-D 

Z - ~~ D~~ 



L> 4 P(Z) 



-<7*(D - - 1)) X + g£> + (2 - q)\u{e' s - 1) 

-(2-^)^(4-1) 
-{2-q)u{e' s -\) 

Da 

X + qD X 

-(7*^(4-1) -(2-9)^(4-1) X + w -1 







-a*u,(4 - 1) 
X 

a*D 



X + u 
2uj 




1 

-1 

X + 2(D-1) 







-IX X + 2D 



= X 



X + qD X 

-{2-q)\u{e' s -l) X + oo -1 

-2X X + 2D 



+ qD 



D 



X 





±07(^-1) X + kj -1 
-IX X + 2D 



= X{(X + g£>)(X + u){X + 2L>) + (2 - g)-cj(4 - 1)X(X + 2D) - 2X(X + g£>)} 



+g J D{D(X + u)(X + 2D) - 2DX - -w(4 - 1)X(X + 2D)} 

X 4 + [2D-2 + ue' s + (D - ^oj(e' s - l))q]X 3 + D[2uje' s + {3D-2 + ^lo 
+D 2 [{2D -2 + Auj- ue' s )q]X + D 3 [2uq] 
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= X A + [2(5 + toe' s ) + (1 + 6+ ^co)q]X 3 + D[2ue' s + (1 + 36 + ^uj)q]X 2 
+D 2 [(26 + 4w)q]X + D 3 [2Loq}. 
The characteristic polynomial is proportional to 

O (Z) = [l + 5+^u£' s ]Z*-[A + 25-(l+5+^Lu)q}Z 3 + [6-ue' s + (u-2)q}Z 2 

_ [ ±-28-(l-6+ 1 -u)q}Z+[l-8+ 1 -ue' s }. 

4.2.3 Von Neumann Analysis 

We successively compute 

MZ) = 26{[2 + uje' s }Z 3 - [6 + ue' s - (2 + ^(e' s + l))q]Z 2 + [6 - ue' s - (2 - u)q]Z 

-[2-^ + ^(4-1)?]}, 
MZ) = 48 2 co{[8e' s -(e' s -l)(2-ue' s )q-^(e' s -l) 2 q 2 ]Z 2 

-[lQe' s -8e' s q + (e' s -l)(l-^)q 2 ]Z 

+ [8e ' s _ (4 _ i)( 6 + U£ ' s)q _ (4 _ if (i _ ±u>)q 2 }}, 

MZ) = 4,5V (4 - 1M4 -g) (2 + u,4) x 

x{[324 - 16(4 " 1)« + (4 " 1)(2 + ^ 2 ]^ - [324 - 1649 + (4 - 1)(2 - ^)g 2 ]}. 

We see that the specific cases q = and 4 = 1 which make 03 vanish shall be treated separately. The 
general case is treated by first checking that |0o(O)| < |0q(O)|, which is obvious. We then notice that 
0i(O) > 0. The relation 0i(O) < 0f(O) is equivalent to — uj(e' s — l)q < 8, which always holds. As for 
relation — 0i(O) < 0*(O), it can be cast as uj(e' s — l)q < 4we' s , which holds true if q < 4. We therefore 
have |0i (0)| < |0*(O)| for q G]0, 4]. We carry on by studying the sign of 

0* 2 (O) = 6 2 u[Ae' s - (4 - l)g][8 + (4 - 1)H > 

and therefore we have to check that 02(0) + 02(0) > and 0^(0) — 02(0) > 

0^(0) + 2 (O) = 2£V(4-l)^ + 32 + 2(4-l)(4-9) 2 ] > 0, 
0^(0) -0 2 (O) = 2«5V(4-l)(2 + ^)(4-9)>0, 

if q G]0, 4[. Last let us study 03 

0^(0) = 45V(4-l)9(4-9)(2 + 4^)[32 + (4-l)(2(4-9)+u,9 2 )]>0, 
0*(O)+0 3 (O) = 8 ( 5V(4-l)9 2 (4-9)(2 + 4^)[8 + ^(4-l)'/] >0, 
0^(0) -0 3 (O) = 86W(e' s -l)q(A-q) 2 (2 + e' s u)[Ae' s -(e' s -l)q]>0. 

Hence we show that 03 and therefore 0o is a Schur polynomial if q G]0,4[ and there remains to treat 
the specific cases. 

4.2.4 Case 9 = 

The case 9 = has once more to be treated separately. We have 

MZ) = (Z-l) 2 [(Z-l) 2 + 6(Z 2 -l) + ^e' s (Z + l) 2 }. 



15 



The corresponding amplification matrix is 



/ 1 
















D-u(e' K -l) 


a; 


-1 




D 


D 


1? 







D 


D-u> 
D 


1 

D 






M4-1) 

D 


-2uj 
D 


2-D 
D 


) 



Anew 1 is a double eigenvalue in two distinct eigensubspaces. The other factor of the characteristic 
polynomial, namely 

MZ) = [1 + 6 + \loe' s }Z 2 - [2 - ue' B ]Z +[1-5 + ±ue' s ], 

should be a Schur (or a simple von Neumann) polynomial. We do have IV'o(O)! < IV'o (^)l an d we 
compute 

^ x {Z)=U{[1 + \us' s ]Z-[1-^e' b ]}. 

Both remaining eigenvalues have a strictly lower to 1 modulus and iterates of the amplification matrix 
are bounded. This holds even if e' s = 1. 

4.2.5 Case q = 4 

In the case when (7 = 4, 

(j) (Z) = [1 + 5 + ^toe' s ]Z 4 + [25 + 2uj]Z 3 + [-2 - coe' s + 4uj]Z 2 + [-25 + 2u]Z +[1-5 + hjje' s ] 
= (Z + 1) 2 {[1 + 5 + \^\Z 2 - 2[1 - U + l ~ue> s ]Z + [1-5 + ±coa' s }}. 
We have a double root Z = —1. We therefore have to study the amplification matrix which reads 



/ 


1 


—a 
















w 


-1 




V 


<r*J«(ei-l) 
D 


D 

-2|«(ei-l) 
D 


D 

D-ui 

D 
-2u 


D 
1 

D 

2-D 


) 


D 


D 


D 


D 



Only the vector (<r, 2,0, 0)* is an eigenvector associated to the eigenvalue —1, and we have increasing 
iterates for G, whatever the study of the other factor of the characteristic polynomial. The value q = 4 
gives rise to instabilities. 

4.2.6 Case e' s = 1 

If e' s = 1, the polynomial 03 is identically zero and we shall study <f>2 for q £]0, 4[ 

fa(Z) = 325 2 lu{Z 2 - [2 - q}Z + 1}. 

This polynomial has two distinct complex conjugate roots with unit modulus. Polynomials <j>2 an d 
therefore 4>q are both simple von Neumann polynomials. 

4.2.7 Synthesis for the Lorentz Kashiwa Model 

The scheme (I14p for the one-dimensional anharmonic Maxwell-Lorentz equations is stable with the 
condition 

5t < 5x/coo. 
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4.3 Young Model 
4.3.1 Model Setting 

The Maxwell system is closed by a discretization of 



namely 



— (P 1 
1 6V > 

1 ( 7" n+ 5 



n+1 



J n +i 



J 



? {J 3 



+ J" 5 )+a; 1 2 ( es -e 00 )eo^-^" 



The explicit version of system 



I)-(I15I) deals once more with the variable 



i+i 



and reads 



3+1 

[1 + 5}£] +1 



B. 



3+1 



81 



[1 + 5\P 



n+1 



[1 + 5 - 2u,(4 - + A 2 [l + 5](£:7 +1 - 28? + fJL x J 
+2w^-[l-<5]j/~', 

[1 + 5- 2oj)V? + 2u(s' s - 1)8] + [1 - 5}J"~K 
[1- 5]^^ +2^(4-1)^-2^, 



3 + 



from which stems the amplification matrix 

/ 1 

G = 



a" 





— CJ 










(l-g)(l+(5)-2wa 


2uj 


1-5 




1+5 


1+5 


1+5 




2uja 


l+5-2u 


1-5 




1+5 


1+5 


1+5 






-2uj 


1-5 


) 


1+5 


1+5 


1+5 



4.3.2 Computation of the Characteristic Polynomial 

The characteristic polynomial of G is 

Z - 1 



P(Z) 



-cr 





a 

(l-q)(l+5)-2Kjq 

1+5 
2tjja 

1+5 
2tjja 

1+5 



z 











2w 


1-5 




1+5 


1+5 r 




l+5-2fj 


1-5 




1+5 


1+5 




2lu 




5 


1+5 


*-T+ 


5 



Therefore setting X = (1 + 6)(Z - 1), 



(1 + 5) 4 P(Z) 





(l + <5)<7 





a*(i + <y) 


X + (1 + S)q + 2uxx -2uj 1 - 


-5 





—2ua 


X + 2to -(1 


-5) 





—2ua 


2uj X + 25 




(1 + S)a 







a*(l + 5) 


X + (l + 5)q 


X 







—2ioa 


X + 2uj -(1-5) 










-X X+l+5 
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X 



X+(1+S)q X 
-2toa X + 2co -(1-6) 
-X X + l + 5 



1 

+ (l + 5) 2 q -2ua X + 2lo -(1 - 5) 
-X X+l+5 

X{X 3 + [25 + 2ue' s + (1 + 5)q]X 2 + [2(1 + 5)coe' s + 2(1 + 5){5 + u)q]X 

+ [2(1 + 5) 2 0Jq]} + (1 + 5) 2 q{X 2 + [2(5 + uj)]X + [2(1 + 5)u]} 

X A + X 3 [25 + 2ue' s + (1 + 5)q] + X 2 [(l + 5)(2ue' s + (1 + 35 + 2u)q)} 

+X[(1 + 5) 2 q(25 + 4a;)] + [2(1 + 5fuq]. 



The characteristic polynomial is proportional to 



[1 + 5}Z A 
-[4-2(5- 



- [4 + 25 - 2uje' s - (1 + 5)q]Z 3 + 2[3 - 2uje' s + (u - l)q]Z 2 
2ue' s -(l-5)q\Z + [l-5]. 



4.3.3 Von Neumann Analysis 

We successively compute 

(f)i(Z) = A5{Z 3 - [3-ioe' s -q}Z 2 + [3-2ioe' s + (uj-l)q}Z- [l-ue' s ]}, 

Mz) = (4<5)M[4(2-^'s)]^ 2 -[24(2-^s) + (4(^-i)-i))'z]^ 

+ [e' s (2-u;e' s )-(e' s -l) q ]}, 
MZ) = (45) 4 qu; 2 (e' s -l){[2e' s (2-u;e' s )-(e' s -l) q ]Z 

-[24(2 - ue' s ) - (4 - l)q - (2 - ue' s )q]}- 

We see that we shall once more treat the cases q = and e' s = 1 separately since 03 is identically 
zero. Let us check the conditions in the general case. First, |</>o(0)| < |<^o(0)l clearly holds as well as 
|0i (0)| < |0i(O)| under the condition u < 2/e' s . If e' s > 1 and q / 0, the condition |<^ 2 (0)| < \M°)\ is 
equivalent to 

(4 -l)q< 24(2-^4). 

If the worst case is q = 2, we must have (4~ 1)2 < 2e' s (2—ue' s ), which is equivalent to oj < (4 + l)/4 2 - 
If the worst case is q = 4, we must have (4 — 1)4 < 24(2 — oje' s ), which is equivalent to u < 2/e' s 2 . 
We wait until the study of 03 to choose between q < 2 and q < 4. The root of 03 is 

_ 24(2 - 1^4) - (4 - l)g - (2 - ue' s )g 

24(2 -^' s ) -(4 -i)g 

The denominator is positive under the same assumption found to ensure |02(O)| < 1 01(0)1. If we want 
|Z| < 1, the condition is hence 

(2 - ue' s )q < 44(2 - uje' s ) - 2(4 - l)q. 

If the worst case is q = 2, we must have 2(2 — we'J < 44(2 — <^4) — 4(4 — I)' w hich is equivalent to 
to < 2/(24 ~ !)• If the worst case is q = 4, we must have 4(2 — ue' s ) < 44(2 — oje' s ) — 8(4 — 1), which 
is equivalent to 4a>4(4 — 1) < 0, which is false. We therefore choose to take q < 2 and the successive 
conditions found are 

2 4 + 1 2 

LO < -7T- , UJ < 



10 < — . 

Co 



-/2 



2e' - 1 
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The more restrictive condition that we have encountered is u < 2/(2e' s — 1), this is therefore our final 
stability condition in addition to q < 2, for which 0o is a Schur polynomial. 

Which are the limiting case we have to study? The three conditions are equivalent if and only if 
£g = 1. In this case, if u has its limit value u = 2 and q its limit value q = 2, we have (p2(Z) = 0. If 
£g 7^ 1, q = 2 and u = 2/{2e' s — 1), then the modulus of the root of 03 is 1, and we conclude that 03 
and therefore 0o are simple von Neumann polynomials. 

4.3.4 Case q = 

If q = 0, the characteristic polynomial has the double eigenvalue 1 

0o(Z) = (Z- 1)\(Z - l) 2 + 5(Z 2 - 1) + 2uje' s Z). 
The corresponding amplification matrix is 



G = 



( 1 
















l+5-2wa 


2uj 


1-5 


1+5 


1+5 


1+5 







2uja 


l+5-2u 


1-5 




1+5 


1+5 


1+5 




V 


2uia 


-2lu 


1-5 
1+5 


) 



The double eigenvalue operates on two distinct eigensubspaces. We have to study the other factor of 
the polynomial 

MZ) = [1 + S]Z 2 - [2(1 - ue' s )]Z + [1-6]. 
We clearly have 1^0(0)1 < IV'oWI- Moreover we compute 

MZ)=4HZ-[1-uje' s }}. 

We recover the condition co < 2/e' B , under which we have a Schur polynomial. For uj = 2/e' s , we have 
a simple von Neumann polynomial (—1 is a root), which allows to conclude. 

4.3.5 Case e' s = 1 

In the case when e' s = 1, only the lasts steps have to be considered. The only problems are the 
condition 1 02 (0) I < 1 02 (0) | and a vanishing ^3. In this specific case, 

<fo{Z) = (45) 2 w(2 - co){Z 2 - [2 - q]Z + 1}. 

If to < 2/e' s , i.e. to < 2, this polynomial is identically zero, and 

<// 2 (Z) = (45) 2 uj(2 - uj){2Z - [2 - ?]}. 

Polynomial (j)' 2 is a Schur one if q G]0, 2] and hence 4>q is a simple von Neumann polynomial. 
If uj = 2, polynomial 02 is identically zero and we compute 

MZ) = 45{Z + \){Z 2 - [2 - q]Z + 1}. 

For q g]0, 2] the roots of Z 2 — [2 — q]Z + 1 are complex conjugate, distinct, and their modulus is 1. 

4.3.6 Synthesis for the Lorentz Young model 

The scheme (|9l)~ (ll5p for the one-dimensional anharmonic Maxwell-Lorentz is stable with the condition 

5x 2 \ 



5t < min 
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5 Harmonic Lorentz Type Media 

Harmonic Lorentz type media are treated thanks to the three above mentioned schemes. The com- 
putation of the amplification matrices and the characteristic polynomials remains unchanged. The 
harmonicity v = is expressed by the parameter 5 = 0. This vanishing value makes 4>\ identically 
zero for the three schemes and the above analysis breaks down. We resume to the analysis of the three 
schemes. 

5.1 Joseph et al. Model 
5.1.1 General Case 

Since 4>\ is identically zero, we want to apply Theorem [2] and study the derivative polynomial of 4>o, 
which we denote by 

ipo(Z) = [4 + 4uje' s }Z 3 - [12 + 6uje' s - 3(1 + co)q]Z 2 + [12 + 4we' B - 4q]Z - [4 + 2ue' s - (1 + u)q]. 

We notice that V'o(O) > an d "00(0) = — (4 — q) —uj(2e' s — q) < for q < 2. We therefore have to check 
that — V'o(O) < V'o(O)) which is equivalent to — (u + l)q < 2loe' s and always holds. We therefore have 
1^0(0)1 < I^Q (0)|. Let us now compute 

ih(Z) = [4^(4 + 3^) +4(1+ u)(2 + UE' s )q- {l+u) 2 q 2 ]Z 2 
+ [-16^4(2 + ue' B ) + 8(wV 8 - 2)q + 4(1 + u;)q 2 ]Z 
+[4^4(4 + ue' B ) + 4(2 - we' B + 6uj + 3u 2 e' B )q - 3(1 + u;) 2 q 2 ]. 

Anew we have 

^i(O) = 4u;4(4 + 3u;e' s ) + (1 + u)q[A(2 + ue' B ) - (1 + u)q] > 0. 

Instead of studying the sign of V'i(O), we will check that — V>i(0) < ipi(0) and ^i(O) < V'i(O)- The 
relation — ^i(O) < ^(O) is equivalent to 

4we' B (2 + ooe' s ) + (4 - q){\ + ufq + Au 2 (e' s - l)q > 0, 

which clearly holds. There remains V'i(O) < V'i(O) which reverts to 

Auj 2 e' s 2 + 4w(e' s - 2 - ue' B ) + (1 + ujfq 2 > 0. 

Cast like this it is not easy to conclude, but we can write it has a polynomial of the variable u 

cu 2 (2e' s - qf + 2^(2(4 - 2) + q) + q 2 > 0. 

The reduced discriminant of this polynomial is 

A' = -8(e' s -l)q 2 (2-q) < 0, 

if < q < 2 and e' s 7^ 1. The polynomial (in uj) is therefore always positive, which we were seeking. 
Hence we have |V , i(0)| < I V^i (0) I an d we can profitably carry on with the computation of ip2(Z) which 
is the product 

ifj 2 (Z) = 8[4w 2 4 2 + (4we' s - Auj 2 e' s - 8to)q + (1 + w)V] x 

x{[4w 2 4 2 + 8coe' s + (4 + 4we' s + 8iu)q - (1 + uj) 2 q 2 ]Z 

-[4cA' s 2 + 8cue' s + (4 - 2ue' B )q - (1 + uj)q 2 )} 
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We notice that 

8[4u;V s 2 + (4we' s - Alo 2 e' s - 8co)q + (1 + tufq 2 } = 8[(4e' s - q)u - q] 2 + Wu;q[4(e' s - 1) + 2e' s ] > 

and we simplify by this factor denoting 

i) 2 {Z) = [4wV s 2 + 8ios' s + (4 + Auj 2 e' s + 8uj)q-{l+Lo) 2 q 2 ]Z 
-[Alo 2 e' 2 + 8ue' s + (4 - 2ue' s )q - (1 + u)q 2 \. 

We see that 

$(0) = 4^ 2 [4 2 + (4 - l)q] + (1 + u) 2 q(4 - q) > 

and to prove |-0 2 (O)| < |$J(0)|, we only have to check that $j(0) + ^ 2 (0) > and $J(0) - ^ 2 (0) > 0. 
We "notice" that 

r 2 (0)+M0) = coq[(4e' s -q)uj + 2(8-q)}>0, 

r 2 (0) - MO) = [2^4 + (1 + w)?][(44 - g)w + 2(4 - g)] > 0, 

which ends the proof in the general case for 5 = 0. 

5.1.2 Case q = 

In the case when q = 0, along with the fact that 1 is a double root "which does not cause any trouble" , 
4>o has the same roots as those of [1 + lue' s ]Z 2 — 2Z + [1 + u>e' s ] which are complex conjugate, distinct, 
and their modulus is 1, if e' s > 1. 

5.1.3 Case q = 2 

The same holds for q = 2 and this time, along with the roots ±i, we have the roots of the polynomial 
[I+oje'^Z 2 — 2{1+uj(e' s — 1)]Z + [l+u;4] which has two complex conjugate, distinct roots with modulus 
1, ife' s >l. 

5.1.4 Case e' s = 1 

Finally if e' s = 1 (and q g]0, 2[), we shall return to the polynomial </>o which can be cast as 

(f> (Z) = [Z 2 - (2 - q)Z + 1][(1 + u)Z 2 - 2Z + (1 + u)]. 

Each of the second degree polynomials has two distinct complex conjugate roots. We therefore have a 
simple von Neumann polynomial except in the particular case when the two polynomials are propor- 
tional and have the same roots, which are then double roots. This is reached if (2 — q) = 2/(1 + lo), 
namely q = 2lo/(1 +uj) or equivalently uj = q/{2 — q). In this case, von Neumann analysis is not useful 
anymore and we have to revert to the amplification matrix 



G = 



( 


1 


—a 





o\ 





-2q + 2 


-l 


q 







1 








V 


a* 


-q 





i / 



The two double eigenvalues of this matrix are (2 — q ± i^J 'q(4 — q))/2 and their each only have one 
corresponding eigenvector 

/ Vg(4-g) g(3-g)=R(2-g)\/g(4-g) -g±yg(4-g) V 
^' 2 2 2 ] ' 
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For each eigenvalue the associated minimal eigensubspace is therefore two-dimensional, which corre- 
sponds to an unstable case. We can say that the scheme is stable for q 6 [0, 2u/(l + uj)[. If we rewrite 
this in physical variables, we have 

2 St 2 . 2 f 2cu 2 St 2 



4cl— ^sin z - < 



00 5x 2 V 2 + 

If 5x 2 < Ac^/uj 2 , this is not a bound on the time step. If 5x if large enough, the bound on the time 
step is 

Sx 2 2 



8t 2 < 



2c^ lo 2 ■ 



5.1.5 Synthesis for the Harmonic Lorentz Joseph et al. Model 

The scheme ([8|)- (jl3p for the one-dimensional harmonic Maxwell-Lorentz equations is stable with the 
condition 



St < ^ if e s > £oo and 5t < J ^ \ if e s = e^, 



this last condition being meaningful only if 5x > 2coo/lui. It is therefore advisable not to use the 
Lorentz- Joseph et al. scheme in the harmonic case for e s = Eoq. 

5.2 Kashiwa Model 
5.2.1 General Case 

Anew the polynomial 4>\ is identically zero and we study the derivative of polynomial 0o, which we 
denote 

MZ) = [4 + 2uje' s ]Z 3 - [12 - 3(1 + ^uj)q]Z 2 + [12 - 2u>e' s - (4 - 2uj)q]Z - [4 - (1 + -u)q]. 

The condition |V>o(0)| < I V^o CO) | is equivalent to (8 — q) + ^u>(Ae' s — q) > which holds for q < 4. Then 
we compute 

MZ) = [4(4 + ^K + 4(2 + .) ? -(l + i w )V]2 2 

-[32we' s + (16 -8u- 8uje' s - 4w 2 £ f s )q + (to 2 - i)q 2 ]Z 
+ [4(4 - uje' s )ue' s + 4(2 - 2oje' s + 5co)q - 3(1 + \u) 2 q 2 ]. 

We check that 

#r(0) = ^[4ue f B + (2 + u)q][16 + 4ue! a -(2 + u;)q] > 0, 
^(0)+Vi(0) = q(A-q)(u; + 2) 2 + Aq(E , s -l)LU 2 + 8E' s (A-q) + 8q>0, 
^(O)-Vi(O) = -q)uj-2q} 2 + 32^(4 - 1)} > 0, 

under the only condition that q G [0,4], which ensures that |V>i(0)| < (^(O)]. Last we compute ip2(Z) 
which can be cast as 

MZ) = i{[(-8 + 44 + q)u + 2q} 2 + 16a(4 - q)}i> 2 {Z) 

with 

MZ) = [(2 + uj) 2 q(A -q)+ 4((4 - l)co 2 q + 2(4 - l)w(4 - g) + 8)]Z - [8 - (2 + w)?]^ + (2 - 
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We check that 



r 2 (o) > o, 



V> 2 *(o) + v> 2 (o) 
$(o)-&(o) 



qu[(6e' B -q)u + 2(8-q)}>Q, 
2(4-q)[4e s u + (2 + u)q] > 0. 



These two last inequalities are strict only if q e]0,4[, and we then have |^2(0)| < 1^2 (0)1- The case 
e' s = 1 is not specific in this general study. 

5.2.2 Case q = 

To treat the specific case when q = 0, we have to revert to the study of 4>o which is here 



The double eigenvalue Z = 1 is the same as in the anharmonic case and does not cause any trouble 
either (minimal eigensubspaces are still one-dimensional). The other factor of the polynomial has 
clearly two distinct complex conjugate roots of modulus 1. This configuration corresponds to a stability 
case for the scheme. 

5.2.3 Case q = 4 

The analysis performed in the anharmonic case remains valid here. The eigenvalue Z = — 1 is double 
and the associated minimal eigensubspace is two-dimensional. Iterates of the amplification matrix are 
therefore linearly increasing and the case is unstable. 

5.2.4 Synthesis for the Harmonic Lorentz Kashiwa Model 

The scheme (|14|) for one-dimensional harmonic Maxwell-Lorentz equations is stable with the condition 



5.3 Young Model 
5.3.1 general Case 

Once more, the polynomial (p\ is identically zero and we study the derivative of polynomial 4>o, which 
we denote 

ipo(Z) = 4Z 3 + [-12 + 6ue' s + 3q]Z 2 + [12 - 8ue' s + 4(w - \)q\Z + [-4 + 2ue' s + q]. 

The condition |?/>o(0)| < IV'oWI ^ s equivalent to (4 — 2oje' s ) + (4 — q) > 0, which we assume (lj < 2/e' s 
and q £ [0, 2]). We carry on computing 

^i(Z) = [{2uje' s + q){8 - {2ujs' s + q))]Z 2 + [{2ue s + q){4{2uoe' s + q) - 16 - Auq) + 16uq]Z 



[1 + -<]Z 4 - 4Z 3 + [6 - coe'^Z 2 -4Z + [1 + -ue' s ] 
(Z - 1) 2 {[1 + \^\Z 2 - 2[1 - ±ue' s ]Z + [1 + ±ue' B )}. 



5t < 5x/c ( 



+ [(2coe' s + q)(8- 3(2ue' B + q)) + Wujq]. 



We see immediately in this formulation that V'i(O) > 0- Besides 



^(0)+Vi(0) 
^(0)-^(0) 



4{-4aA' s 2 + 8^e' s + 4(cj(e' B + 1) + l)q - q 2 }, 
2{[2co - qf + 4uj(e' s -!)[? + u>& + 1)}] > 0, 
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Let us note f(u>) = —Auj 2 e' 2 + 8uj£' s + 4(cj(e' s + 1) + l)q — q 2 , we want to prove that this quantity is 
positive for uj g]0, 2/(24 — 1)]. We derive to obtain f'(oj) = —4:e' s q — 8uje' s 2 + 8e' s + Aq which vanishes 
at u = (2e' s — (e' s — l)q)/2e' s 2 . This corresponds to values of u between l/e' s 2 (value for q = 2) and 
1/4 (value for q = 0), which always belong to the interval ]0, 2/(2e' s — 1)]. At this point we have a 
maximum of the function f(u>). To conclude, we only have to evaluate the limit for f(ui) as uj — > 
and the value of f(uo) at uj = 2/(2e' s — 1). If both values are positive, f(oj) will be positive on the 
whole interval. 

lim/(w) = 4q-q 2 >0, 
7(2/(24 - 1)) = (24 - 1) 2 (2 - g) + 2(24 -1)+ 4(4-1) (2e' s + 1)>0, 

If 4 = 1) q = 2 and uj = 2, we have V'i(O) = — "01 (0)- We will treat this case apart. We finally compute 
ip2{Z) which can be cast as 

MZ) = 8[(2coe' s - q) 2 + 8u;(4 - l)q]fc(Z) 

with 

MZ) = [-4w 2 4 2 + 8ue' s + 4(w(4 + 1) + l)g - q 2 ]Z - [{2ue' s + (1 - w)g)(-4 + 2ue' s + g)]. 

We notice that ip2 is identically zero if e' s = 1 and g = 2u;e^ = 2oj, which has to be treated separately. 
In the opposite case, we check that 

^l(o) > /H>o, 

^2(0) +^ 2 (0) = qu[8 - 2coe' s - q] > 0, 
V4(0)-M0) = [2a;4 + g][8-4a;4 + (a;-2)g]. 

The quantity 8 — 4cj4 + (w — 2)q is minimum if q = 2 (since to < 2) and is then equal to 4 — 2u;(24 — 1)- 
As in the anharmonic case if uj < 2/(24 — 1) this quantity is positive, and if u = 2/{2e' s — 1) this 
quantity is zero. Yet we want to show that -00 is a Schur or a simple von Neumann polynomial, we 
therefore have to revert to the study of 4>q. Once more the cases e' s = 1, q = 2 and u> = 2 have to be 
treated specifically. 



5.3.2 Case q = 

In the case when q = 0, we revert to the direct study of 0o 

<f) (Z) = Z 4 - [4 - 2uje' s ]Z 3 + 2 [3 - 2ue' s ]Z 2 - [4 - 2ue' s ]Z + 1 
= (Z-1) 2 [Z 2 -2(1-„4)Z+1], 

which has Z = 1 as a double root, which is no more a problem as in the anharmonic case. Both 
other roots are complex conjugate, distinct and have a unit modulus if u>e' s < 2. If e' s = 1 and uj = 2, 
Z = — 1 is also a double root. Then we have 



G + Id 



/ 2 

a* 


V 





1 
1 



\ 
_ 1 

I 2 

1 

1 / 



which has only (a, 2, 0, 0)* as eigenvalue (associated to 0). This is a cause of instability for the scheme. 
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5.3.3 Case q = 2 

The case q = 2 is treated by the general case except when u> = 2/(2e' s — 1). In this case 2uje' s = 2 + u, 
hence 

00 (Z) = Z 4 + loZ 3 + 2[uj -1}Z 2 +uoZ + 1 
= (Z + l) 2 (Z 2 -(2-co)Z + l). 

We therefore have to study the stable subspaces through the amplification matrix for the eigenvalue 
— 1. We have 

( 2 -a \ 

a* uj-2 2tu -1 
2-w 2-2w 1 
\ 2-oj -2uj 2 ) 

which has a unique eigenvector (associated to 0), namely (a, 2, —1, —2)'. This is an unstable case. 

5.3.4 Case e' s = l 

Only the case when e s = 1, q = 2u> remains to study, in which case ip2 vanishes. We compute 4>q which 
is equal to 

MZ) = (Z 2 - 2[1 - uj]Z + l) 2 . 



G + Id 



The two complex conjugate roots 1 — uj ± iy / uj(2 — uj) are both double with modulus 1. We there- 
fore have to study the associated stable subspaces. The only associated eigenvectors are (a, u> =F 
i-\/uj{2 — uj)) 1 respectively and the associated minimal eigensubspaces are two-dimensional, which leads 
to instability. If e' s = 1, we should assume q < 2u, which in physical variables reads 5x > 2coo/wi 
which is not a stability condition. It should therefore be avoided to use the Lorentz- Young scheme in 
the harmonic case for e s = when q can reach the value 2u, i.e. if uj > 1. Another way to see this 
condition is to give u < 1 as a stability condition if e' s = 1. 

5.3.5 Synthesis for the Harmonic Lorentz Young et al. Model 

The scheme ([9| M[15p for the one-dimensional harmonic Maxwell-Lorentz equations is stable with the 
condition 

. ( Sx 2 \ . r , . 

ot < mm — = — , , if e s > e^, and ot < mm 




\Z2coo' ujiy/2e' s - 1 

6 Basic Polynomials in Dimension 1 

The previous computations lead us to define basic polynomials associated to each one-dimensional 
scheme. We will see that these polynomials will prove useful in higher dimensions. 

Debye (Joseph et al.) 

Pdj(Z) = [1 + 5e' s ]Z 3 - [3 + Se' s - (1 + 5)q]Z 2 + [3 - 6e' s - (1 - b)q\Z - [1 - fe'J 
= [1 + 5e' s ]Y 3 + [2Se' s + (1 + 5)q]Y 2 + [(1 + 35)q}Y + [2Sq] . 

Debye (Young) 

P DY (Z) = [(1 + 5a)(l + 5)]Z 3 -{3 + 5 + 5a + 35 2 a-(l + 5)q]Z 2 
+ [3 — 6 — Sa + 35 2 a - (1 - S)q]Z - [(1 - 5a)(l - 5)} 
= [(1 +S)(1 + Sa)]Y 3 + [25(1 + a) + (1 + 5)q]Y 2 + [(1 + 35)q]Y + [2Sq\. 
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Lorentz (Joseph et al.) 



P LJ (Z) = [1 + 5 + uje' s \Z 4 -\A + 25 + 2ue' s - (l + 5 + u)q]Z 3 + [6 + 2coe' s - 2q]Z 2 
-[4 - 25 + 2ue' s -(1-5 + uj)q\Z +[1-5 + gje' s ] 
= [1 + 5 + coe' s ]Y 4 + [25 + 2ue' s + (1 + 5 + u;)q}Y 3 + [2ue' s + (1 + 35 + 3co)q]Y 2 
+ [2(5 + 2uj)q}Y + [2uq\. 



Lorentz (Young) 

Ply(Z) = [1 + 5}Z 4 - [A + 25-2uje' s -(l + 5)q\Z 3 + 2[3-2uje' s + (uj-l)q\Z 2 
-[4 - 25 - 2toe' s - (1 - 5)q]Z +[1-5] 
= [1 + 5}Y 4 + [25 + 2u;e' s + (1 + 5)q]Y 3 + [2ue' s + (1 + 35 + 2uj)q]Y 2 
+ [2(5 + 2uj)q}Y + [2uq\. 

7 The Two-Dimensional Space Case 
7.1 Maxwell Equations 

In dimension 2, the field may be decoupled into two polarisations which lead to different schemes 
and also a different number of variables. We use here similar notations as those introduced in the 
one-dimensional case, namely X x = Coo5t/5x, X y = CooSt/Sy, a x = X x (e l ^ x — 1), a y = X x (e l ^ v — 1), 
q x = Wx\ 2 and q y = \(J y \ 2 . 

7.1.1 Polarisation (B x , B y , E z ) 

The polarisation (B x , B y , E z ) is also called the transverse electric polarisation TE Z . 



Lorentz (Kashiwa et al.) 



Plk(Z) 





A y \ c z,j,k+l c z,j,k) > 
A x ( £ z,j+l,k - £ z,j,k) > 



u z,j,k u z,j,k 



x x (V + ^ -s n+ ^ 



) 
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7.1.2 Polarisation (B z ,E x ,E y ) 

The polarisation (B z , E x , E y ) is also called the transverse magnetic polarisation TM Z . 



B 



' + 2 _ R n 2 — _\ (cn _ cn \ , \ fen _ cn \ 



2 ' ^ 2 



n+ 



2 ■" i 2 z '-?^"2 2 



2J + ifc + 



7.2 The Debye— Joseph et al. Scheme 
7.2.1 Polarisation (B x ,B y ,E z ) 

Coupling polarisation (B x , B y , E z ) with the Debye-Joseph et al. scheme, we obtain the amplification 
matrix 



G = 



( 



\ 



1 



1+5^ 





1 

1+54 
— cr,. 



On 



(l-8e' B )-(l+S)(q x +q y ) 

1+54 

-(<?x + %) 
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1+54 
1 



\ 



associated to the variable (B 2 1 ,B 



cn 

-i c 'z. 



z ■ k , T>™ ■ k y. The characteristic polynomial is propor- 



tional to the characteristic polynomial in dimension 1 for the same scheme with 1 as an extra root 

MZ)=YP DJ (Z). 

In polynomial Pdj(Z), the variable q means q = q x + q y in the two-dimensional case. The polynomial 
only depends on this sum and not on the separate values of q x and q y . 

The general case treated in dimension one concludes to a Schur polynomial, we therefore have a 
von Neumann polynomial here. We also see easily on the amplification matrix in the case g = (if 
q x only or q y only vanish, we do not have a specific case), that the eigensubspaces associated to the 
eigenvalue 1 are indeed stable. As for the particular case e' s = 1, which gave rise to two complex 
conjugate eigenvalues, different from 1, we may add this new eigenvalue with the same conclusion, 
namely stability with the condition q = q x + q y < 4 i.e. \^2coo5t < 5x if 5x = 5y. 



7.2.2 Polarisation (B z ,E x ,E y ) 

Coupling polarisation (B z , E x , E y ) with the Debye-Joseph et al. scheme yields the amplification 
matrix 



G 



V 



1 

1+54 

0±gK 
1+54 



(1-54)- (1+5)^ 
1+54 

-Qy 
1+54 





25 
1+54 

1 






—fx 
(1+5)0^ <J* 

1+54 

<?X<?y 

(l-8e' s )-(l+5)g x 
1+54 

-Qx 






\ 












25 




1+54 




1 


/ 



associated to the variable (B . 2 1 , 1 , £ n 



,j+i,fe+i' x,j+±,k> x,j+±,k> y,j,k+^ y,j,k+± 



cr 



)*. Anew we have a propor- 



tional polynomial to that of the one-dimensional case with two extra roots 

MZ) = Y[(l + Se' s )Y + 26e' s ]P DJ (Z), 
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which are equal to 1 and (1 — 5e' s )/(l + Se' s ) respectively and which are not roots in dimension 1 in 
the general case. 

Only the root 1 could induce stability problems and we have seen that only the case q = would 
make it a multiple root. The matrix G is then block-diagonal, with two rank-two blocks identical to 
that of the one-dimensional case. The stability is once more ensured with the condition q = q x +q y < 4. 

7.3 The Debye— Young Scheme 
7.3.1 Polarisation (B x ,B y ,E z ) 

Coupling polarisation (B x , B y , E z ) with Debye- Young scheme, we obtain the amplification matrix 



G 



( 



\ 



1 



°l 

l+5a 







1 

a* 

l+5a 





-a y \ 

a x 

(l+5)(l-5a)+46 2 a-{l+5)q 1-8 25 

(l+5)(l+<5a) 1+5 l+Sa 

25a 1-6 I 

1+5 1+5 / 



associated to the variable (B ? ,,B 2 i , , £ n ■ , , V n ■ u )*■ The computation of the characteristic 

v x,j,fe+5' y,j + ^,k } 2J,«' z,3,k' * 

polynomial leads to the same polynomial as in one dimension but with 1 as an extra eigenvalue 

MZ) = YPdy(Z). 

Anew we revert to the arguments of the one-dimensional case, the double eigenvalue 1 of the q = 
case not being a problem. 



7.3.2 Polarisation (B z ,E x ,E y ) 

Coupling polarisation (B z , E x , E y ) with the Debye- Young scheme, we obtain the amplification matrix 



G = 



\ 



1 

_* 

y 

l+5a 



l+5a 





Uy 

(l+5)(l-5a)+45 2 a-(l+5)q y 
(l+5)(l+5a) 
25a 
1+5 

l+5a 







1-5 25 
1+5 l+5a 

1-5 

1+5 






CFxCTy 

T+5a 



{l+5){l-5a)+A5 2 a-(l+5)q x 
(l+6)(l+5a) 
26a 
1+6 







1-6 25 
1+5 l+Sa 

1-5 

1+6 



\ 



J 



z,3- 



. L i,f n ..i ,,V n i x i,r., ± i )*. The computation of the 

■h k +V x J+i< k x,3+i,W y,j,k+±' y,3,k+i' * 



associated to the variable (B 

characteristic polynomial leads to the same polynomial as in dimension 1 but with two extra eigenvalues 

MZ) = Y[(l + <5)(1 + Sa)Y + 25(1 + a)]P DY (Z) 

which are 1 and (1 — <5)(1 — 5a)/(l + 5)(l + 5a). Anew we revert to the argument in the one-dimensional 
case, the triple eigenvalue 1 of case q = not being a problem. 



7.4 The Lorentz-Joseph et al. Scheme 
7.4.1 Polarisation (B x , B y , E z ) 

Coupling polarisation (B x , B y ,E z ) with the Lorentz-Joseph et al. scheme, we obtain the amplification 
matrix 



G = 



( 


1 





-Oy 





o \ 







1 












26a* y 


25a* 


2-(l+5+u)q 


l-5+u)e' s 


2u) 




l+5+u>e' e 


l+5+u)e' s 


(l+S+u>e' B ) 


l+5+uje' s 


l + 6+L0£' s 










1 








V 


a y 




-q 





1 / 
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associated to the variable (B . ? 1 ,B . 2 , , , £™ • £™ ■ i, X>™ • .)*. The computation of the character- 

v x,i,k+i y,i + -£,k z ,h K z,3, K 



Mi -J - n „ y 

istic polynomial leads to the same polynomial as in dimension 1 but with the extra eigenvalue 1 

<h(Z) = YP LJ {Z). 
The triple eigenvalue in the q = case is not a problem. 



7.4.2 Polarisation (B z ,E x ,E y ) 

Coupling polarisation (B z , E x , E y ) with the Lorentz-Joseph et al. scheme, we obtain the amplification 
matrix 



G 



( 



1 

l+8+ue' B 


~ a v 

25a* 
1+5+ujs' s 





Oy 

-(l+8+uj)q y 
1 









l-S+uje' B 
l+8+u)e' s 











2uj 


1 









2-(l+5+u)q x 
1 

-Qx 








l-S+we' B 
l+8+u>e' s l+5+uje s 


1 








2u) 



\ 



associated to the variable (£> 



cr, 

1 > C „ 



cn—1 -j-jr, 



cn c 

' C „, „' I. . I ! C , 



n— 1 



The 



characteristic polynomial is once more proportional to the one-dimensional polynomial 

= y[(i + <5 + ^' s )y 2 + 2(5 + ^)y + 2<]p LJ (z) 

= Y[(l + 5 + uoe' s )Z 2 -2Z + (1-S + ue' s )}P LJ (Z). 
In the anharmonic case, and by the von Neumann technique, we check easily that 

^O(Z) = [1 + 5 + u;e' s ]Z 2 - 2Z + [1 - <5 + w^] 
is a Schur polynomial. Besides the double root 1 if q = is still no problem. 

In the harmonic case, tpo(Z) has two distinct complex conjugate roots with modulus 1. This is not 
a problem in itself, except if e' s = 1 in which case ^q{Z) is also a factor in Plj(Z) 

MZ) = (Z-l)[(l + ooe' s )Z 2 -2Z + (l+uje' s )}P L j(Z) 

= (Z - 1)[(1 + ooe' s )Z 2 -2Z + (1+ uje' s )] 2 [Z 2 - (2 - q)Z + 1]. 

This is already the case when we have detected double eigenvalues in dimension 1, giving rise to 
instabilities for q = 2w/(l + u). It is better to avoid this scheme in the case when e' s . 



7.5 The Lorentz-Kashiwa et al. Scheme 
7.5.1 Polarisation (B x ,B y ,E z ) 

Coupling polarisation (B x ,B y ,E z ) with the Lorentz-Kashiwa et al. scheme, we obtain the amplifica- 
tion matrix 
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associated to the variable (B . J , , B . 2 , , , £ " ■ ■ j.)*. The calculation of the characteristic 

V + ^ y,] + ^,k 2 >Ji K ' Z:Ji K ' 

polynomial leads to the same polynomial as in dimension 1 with the extra root 1 

Mz) = YP LK {Z). 
The triple eigenvalue 1 of case q = is not a problem. 



7.5.2 Polarisation (B Z ,E X ,E, 



Coupling polarisation (B z , E x , E y ) with the Lorentz-Kashiwa et al. scheme, we obtain the amplifica- 
tion matrix 
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associated to the variable (B . 2 1 , 1 , £ n . , i , , V n 



The 



'jn cn Tin rjn Xt 

x,j+±,k>°x,j+±,k^y,j,k+^ yj t k+^ U y,j,k+±> ■ 

computation of the characteristic polynomial leads to a polynomial proportional to the one-dimensional 
one 



Mz) 



Y[(l + 5 + -ue' s )Y 2 + (2(5 + coe' s ))Y + (2ue' s )]P LK (Z) 



1 



1 



= Y[(l + 5 + -ue' s )Z 2 - (2 - ue' s )Z + (1-5 + -coe' s )]P LK (Z). 
In the anharmonic case, and by the von Neumann technique, we check easily that 



1 



1 



MZ) = [l + 5 + -ue' s ]Z 2 - [2 - ue' s ]Z +[1-5 + -ue' s ] 

is a Schur polynomial. Besides, the root 1, which is a double one if q = 0, does not lead to any problem. 

In the harmonic case, we have the extra roots 1 and two complex conjugate roots of modulus 1, 
which are not roots of Plk(Z). The stability is hence given under the same conditions as in the 
one-dimensional case. 



7.6 The Lorentz Young Scheme 
7.6.1 Polarisation (B x ,B y ,E z ) 

Coupling polarisation (B x , B y , E z ) with the Lorentz- Young scheme, we obtain the amplification matrix 
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associated to the variable (B . ? , , B . 2 1 , , £ " ■ V™- J71 „■ t 2 )*. The computation of the character- 

v X,J,K+5 1/J + 5,K 2>J) K ' 2J,K ' 

istic polynomial leads to the same polynomial as in dimension 1 but with 1 as an extra eigenvalue 

MZ) = YP LY (Z). 
The triple eigenvalue 1 of the case q = is not a problem. 



7.6.2 Polarisation (B z ,E x ,E y ) 

Coupling polarisation (B z , E x , E y ) with the Lorentz- Young scheme, we obtain the amplification matrix 
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associated to the variable (B 



V 



i • 



V 



The 



.,, , . , , -x,j+i,fc' ' x,3+±,W U x,j+±,k^y,3,k+^ ' y,j,k+^ u y,j,k+^ ' 
computation of the characteristic polynomial leads to a polynomial proportional to that of dimen- 
sion 1 



MZ) = Y[(l + 6)Y 2 + (2(6 + u J s' s ))Y + (2Lue' s )]P LY (Z) 
= Y[(l + 5)Z 2 - (2 - toe' s )Z + (1 - 5)}P LY (Z). 

In the anharmonic case, and by von Neumann technique, we check easily that 

MZ) = [1 + S]Z 2 - [2 - 2ioe' s ]Z + [1 - 5} 

is a Schur polynomial. Besides, the root 1 which is a double one if q = is not a problem. 

In the harmonic case, we have the extra roots 1 and two complex conjugate roots of modulus 1, 
which are not roots of Ply(Z). The stability is therefore ensured under the same conditions as in the 
one-dimensional case. 



8 Conclusion 

We have studied the stability of numerical schemes for Maxwell-Debye and Maxwell-Lorentz equations 
in space dimension 1 and 2. In dimension 2, the characteristic polynomials of each scheme and in both 
polarisation happen to be proportional to the characteristic polynomials for the same scheme in space 
dimension 1. In all the cases, the extension to dimension 2 goes with an extra root 1 compared to the 
one-dimensional case. This is the only extra root in the TE Z polarisation. For the TM Z polarisation, 
there is one other extra root for the Debye equation and two other extra roots for the Lorentz equation, 
all these roots being on the unit circle. For the Yee scheme applied to the raw Maxwell equations, 
the stability condition is q < 4 in dimensions 1, 2 et 3, recalling that q = q x + q y in dimension 2 
(q = max(q x + q y ,q x + q z ,q y + q z ) in dimension 3). The results are gathered in two tables according 
to e s = £oo or not. 

For each model, we have at least one scheme for which the stability condition is the same as for the 
raw Maxwell equations (q < 4). In Young models, the extra conditions correspond to a fine enough 
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Model 


Scheme 




dimension 1 


dimension 2 (&c = Sy) 


Debye 


Joseph et al. 


q < 4 


£i < i£ 

Coo 


St< * x 


Debye 


Young 


q < 4, 6 < 1 


5i < min(^L,2t r ) 

v Coo 7 


bt < min( ** , 2t r ) 


Lorentz 


Joseph et al. 


q < 2 


5i < l x 

V2Coo 


<5t<#- 


Lorentz 


Kashiwa et al. 


q < 4 


<ft < ^ 

Coo 


st< 


Lorentz 


Young 


q < 2, w < ^ 


St < min( , 2 ) 


tfi < min( 9 fa , 2 ) 


Harm. 


Joseph et al. 


9 < 2 






Harm. 


Kashiwa et al. 


g < 4 


Coo 


st< -^r- 


Harm. 


Young 


? < 2, w < 2eJ 2 _ 1 
or 

g < 2, w < 2£ , 2 _ 1 


5i < min( &x , 2_) 


& < min( ^ , I—) 



Table 1: Stability of schemes for e s > e, 



discretization of Debye and Lorentz equations respectively... because stability is not the only issue. 
Applications to classical materials show in general that the condition due to the Maxwell equations is 
the more restrictive one and not conditions due to the constitutive law of the material. 

Computations in dimension 3 are too tedious to be carried out by hand. They have been automated 
(see 0). 
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